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ABSTRACT 

We investigate the hydrodynamics of the interaction of two supersonic winds in binary sys- 
tems. The collision of the winds creates two shocks separated by a contact discontinuity. The 
overall structure depends on the momentum flux ratio r/ of the winds. We use the code RAM- 
SES with adaptive mesh refinement to study the shock structure up to smaller values of r/, 
higher spatial resolution and greater spatial scales than have been previously achieved. 2D 
and 3D simulations, neglecting orbital motion, are compared to widely-used analytic results 
and their applicability is discussed. In the adiabatic limit, velocity shear at the contact discon- 
tinuity triggers the Kelvin-Helmholtz instability. We quantify the amplitude of the resulting 
fluctuations and find that they can be significant even with a modest initial shear. Using an 
isothermal equation of state leads to the development of thin shell instabilities. The initial 
evolution and growth rates enables us to formally identify the non-linear thin shell instabil- 
ity (NTSI) close to the binary axis. Some analogue of the transverse acceleration instability 
is present further away. The NTSI produces large amplitude fluctuations and dominates the 
long-term behaviour. We point out the computational cost of properly following these instabil- 
ities. Our study provides a basic framework to which the results of more complex simulations, 
including additional physical effects, can be compared. 

Key words: hydrodynamics — instabilities — binaries: general — stars: massive — stars: 
winds, outflows — methods numerical 



1 INTRODUCTION 

The stellar winds of massive stars are driven by radiation pres- 
sure to highly supersonic terminal velocities Voo ~ 1000 — 
3000 km s _1 , with mass loss rates that can reach M « 10~ 6 
M yr _1 in O stars and 10~ 4 M yr _1 in Wolf-Rayet stars ( |Puls| 
et al. 2008 ). The interaction of two such stellar winds in a bi- 
nary system creates a double shock structure where the material 
is condensed, heated and mixed with important observational con- 
sequences (see Pittard et al. 2005 for a review). For instance, these 
colliding wind binaries (CWB) have much larger X-ray luminosi- 
ties than seen in isolated massive stars due to the additional emis- 
sion from the shock-heated material. The increased density in the 
shock region also has an impact on the absorption of light within 
the binary. Further away from the system, free-free emission is de- 
tected in the radio, possibly supplemented by synchrotron radiation 
from electrons accelerated at the shock. High-resolution imaging in 
infrared (Tuthill et al. 1999 ) and radio ( Dougherty et al. 2003] has 
made it possible to trace the large scale spiral structure created by 
the winds with the orbital motion of the stars. The interpretation of 
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these observations requires knowledge of the shock structure and 
geometry. 

Assuming a purely hydrodynamical description, the interac- 
tion results in the formation of two shocks separated by a contact 
discontinuity. In the adiabatic limit, the gas behind the shock is 
heated to temperatures T ~ M 2 T W (where T w is the wind temper- 
ature and A4 > 1 is the Mach number of the wind). The structure is 
shaped primarily by the momentum flux ratio of the winds ( |Lebe-| 
|dev&Myasnikov|1990| > 

M2V002 

n = — • (i) 

Mi Voo 1 

The subscript 1 stands for the star with the stronger wind, the sub- 
script 2 for the star with the weaker wind. For reasons of symmetry, 
the contact discontinuity is on the midplane between the stars when 
77 = 1. Pilyugin & Usov (2007) obtained a complete semi-analytic 
description of the interaction region for this specific case. When 
77 7^ 1 the shock structure bends towards one of the stars as the 
stronger wind gradually overwhelms the weaker wind. This leads 
to a bow shock shape close to the binary and the contact discontinu- 
ity shows an asymptotic opening angle at large scales (neglecting 
orbital motion, Girard & Willson 1987). The shock structure must 
then be derived from numerical simulations ( |Luo et aL]|1990) . It 
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depends on other parameters (Mach number, velocity ratio of the 
winds) and, crucially, on the cooling properties of the gas. Cool- 
ing becomes efficient when the radiative time scale of the shocked 
flow b ecomes shorter than its dynamical time scale ( |Stevens et al.| 
1992). In this case, the kinetic energy of the wind (typically ~ 10 
erg s _1 ) is radiated away and the incoming gas is strongly deceler- 
ated at the shock (v = Voo / M 2 in the isothermal limit compared to 
v = Voo/ 4 in the adiabatic limit). The interaction region becomes 
thin and the double shock structure indistinguishable from the con- 
tact discontinuity. Analytical solutions for the interaction geome- 
try can be derived in the limit of an infinitely thin shock ( Girard & 
Wills on|1987||Luo et al.|1990||Dyson et alTTT993| [Canto et al.|1996[ 



2 NUMERICAL SIMULATIONS 



Gayley 2009 see §3 below). 



The analytical solutions provide useful approximations but 
their validity may be questioned as numerical simulations show 
that shocks become unstable (see §4). The contact discontinuity 
separates two media with different tangential velocities, triggering 
the Kelvin-Helmholtz instability (KHI) in adiabatic or radiatively- 



inefficient shocks. The impact is more or less pronounced ( Stevens 



et al.|| 19921 [Lemaster et al.||2007] jPittard 2009, Parkin & Pittard 



2010, van Marie et al. 201 1 ) and has not been quantified yet. Thin 



shocks become violently unstable and have garnered more atten- 
tion. The instability was initially seen in simulations where the 
gas was assumed to be isothermal, mimicking the effect of effi- 
cient cooling ([Stevens et al.|19 92 , Blondin & Koerwer 1998 but see 
Myasniko v et al.|1998| ), and has since also been seen in simulations 
including a more realistic treatment of radiative cooling ( Pittard 
2009, van Marie et al. 2011). The resulting mixing and variability 
can have important observational consequences. The origin of the 
instability remains unclear ( Walder & Folini 1998}. Two mecha- 
nisms have been proposed in the thin shell limit: the non-linear thin 
shell instability (NTSI, Vishniac 1994) and the transverse accelera- 



tion instability (TAI, Dgani et al. 1993 1996) ; both may be at work 
in colliding winds (Blondin & Koerwer 1998). 

Much progress has been made in including more realistic 
physics in simulations of CWB (wind acceleration, gravity from 
the stars, radiative inhibition, cooling functions, heat conduction, 
orbital motion etc.). These are undoubtedly important effects to 
consider when comparing with observations but they complicate 
the comparison with basic analytical expectations which, in turn, 
makes it more difficult to assess their contributions. Here, we 
present simulations neglecting all these effects, assuming a poly- 
tropic gas P oc p 1 with 7 = 5/3 (adiabatic) or 7 = 1 (isothermal). 
Our purpose is to understand how the shock region compares to 
expectations and to constrain the conditions giving rise to instabil- 
ities particularly in the limit of low 77. We performed a systematic 
set of 2D and 3D numerical simulations using the hydrodynamical 
code RAMSES (Teyssier 2002) with adaptive mesh refinement, al- 
lowing us to reach the high resolutions required for thin shocks and 
low 77 while keeping a wide simulation domain to study the asymp- 
totic behaviour (§2). Notable features of the wind interaction re- 
gion are discussed and compared to the analytical solutions: shock 
location, width, opening angle and the presence of reconfinement 
shocks at low 77 (§3). We present our investigations of the instabili- 
ties in the adiabatic and isothermal case in §4. We find that the non 
linear thin shell instability (NTSI) is the dominant mechanism for 
isothermal winds. We then replace this work in its larger context, 
discussing the impact that including additional physics would have 
on our conclusions and the computational cost required to follow 
the instabilities (§5). 



We use the hydrodynamical code RAMSES (Teyssier 2002) to per- 
form our simulations. This code uses a second order Godunov 
method to solve the equations of hydrodynamics 
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where p is the density, v the velocity and P the pressure of the gas. 
The total energy density E is given by 



(5) 



7 is the adiabatic index, its value is 5/3 for adiabatic gases and 
1 for isothermal gases. For numerical reasons 7 is set to 1.01 for 
isothermal simulations (Truel ove et al.|19 98). We use the MinMod 
slope limiter. We compare our simulations with analytic solutions 
in §3. In order to do this, we prevent the development of instabilities 
in the shocked region by using the local Lax-Friedrich Riemann 
solver, which is more diffusive. An exact Riemann solver is used to 
study the development of instabilities in ^4] We perform 2D and 3D 
simulations on a Cartesian grid with outflow boundary conditions. 
We use adaptive mesh refinement (AMR) which enables to locally 
increase the spatial resolution according to the properties of the 
flow. In 2D the grid is defined by a coarse resolution n x — 128 with 
up to 6 levels of refinement. In 3D the grid is defined by n x — 32 
with up to 5 levels of refinement. The refinement criterion is based 
on density gradients. 



2.1 Model for the winds 

Our method to implement the winds is similar to the one developed 
by |Lemaster et a l. ( 2007) and described in the appendix of their pa- 
per. The main aspects are recalled here for completeness. Around 
each star, we create a wind by imposing a given density, pressure 
and velocity profile in a spherical zone called mask. The masks are 
reset to their initial values at all time steps to create steady winds. 
The velocity is purely radial and set to the terminal velocity of 
the wind in the whole mask. Setting the velocity to supposes the 
winds have reached their terminal velocity at the interaction zone. 
This might not be applicable for very close binaries or if 77 <C 1 
because the shocks are then very close to one of the stars. Our 2D 
setup differs from those usually found in the literature (e.g. Stevens 
|et al.|1992[|Brighenti & D'Ercolel^j^lPittard et al.|2006) in that 
we work in the cylindrical (r, 6) plane instead of the (r, z) plane. 
A drawback of our 2D method is that the structure of the colliding 
wind binary is not identical when going from a 2D to 3D simula- 
tion with the same wind parameters. However, as described later, 
we found that the 3D structure is mostly recovered in 2D by us- 
ing the scaling ^ffjsD — V 772D. An advantage of our 2D approach 
is that it is straightforward to include binary motion without re- 
sorting to full 3D simulations. Such simulations will be described 
elsewhere (see Lamberts et al. 2011 for preliminary calculations). 
The density profile is determined by mass conservation through the 
mask 
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where r is the distance to the centre of the mask. The pressure is 
determined using Pp~ J — K with K constant in each region. 
Time is expressed in years and mass loss rates are expressed in 
10 _8 M© yr _1 . We decide to scale all distances to the binary sepa- 
ration a. This way the results of a simulation can easily be rescaled 
to systems with a different separation. For each simulation, the in- 
put parameters are the mass loss rate, terminal velocity and Mach 
number A4 at r = a of each wind. We then derive the hydrodynam- 
ical variables at a. After that the corresponding density, pressure 
and velocity profile in the mask are computed. 

For r\ <C 1 the shocks form very close to the second star. In 
this case, the mask of the star has to be as small as possible so that 
the shocks can form properly (Pitta rd|1998| ). However a minimum 
length of 8 computational cells per direction is needed to obtain 
spherical symmetry of the winds. We thus fix the size of the masks 
to 8 computational cells in each direction for the highest value of 
refinement. We performed tests with a single star for different sizes 
of the mask ranging from 0.03a to 1.5a. The tests were performed 
for n x = 128 and 4 levels of refinement. The resulting density pro- 
files all agree with the analytic solution with less than 1 % offset. 
The surrounding medium is filled with a density pamb — 10 _4 p(a) 
and pressure P am b — 0. LP(a) . This initial medium is pushed away 
by the winds. Simulations with different pamb and P a mb show the 
same final result, to round-off precision. The size of the computa- 
tional domain varies between lb ox = 2a and hox — 80a accord- 
ing to the purpose of the simulation. Except where stated other- 
wise, we took Mi = M 2 = 10~ 7 M yr" 1 , Mi = M 2 = 30, 
Voo2 = 2000 km s _1 and 77 was varied by changing Vooi- Hence, 
our low momentum flux ratios can imply very high velocities for 
the first wind. 
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Figure 1. Density map of the interaction zone for 77 = 1/32 = 0.03125 
(3D simulation). It is a cut perpendicular to the line of centres taken from a 
3D simulation. A zoom on the binary system is shown at the bottom right 
corner. The stars are positioned at the intersections of the dotted lines. The 
first star has coordinates (0,0), the second one has coordinates (a,0). There 
are three density jumps (for increasing x). The first shock separates the 
unshocked wind from the first star from the shocked wind. The contact dis- 
continuity separates the shocked winds from both stars. It intersects the line 
of centres at the standoff point R s . The second shock separates the shocked 
and unshocked parts of the wind from the second star. As the wind from 
the second star is collimated, there is a reconfinement shock along the line 
of centres. R(0±) is the distance between the contact discontinuity and the 
first star, #1 is the polar angle. The asymptotic opening angle is given by 
0ioo. I is the distance to R s along the contact discontinuity. 



3 THE SHOCK REGION 



In this section we study the dependence on 77 of the geometry of the 
interaction zone. We discuss he analytic solutions for the colliding 
wind geometry, in 2D and 3D, to which we compare our simu- 
lations. Simulations are performed with adiabatic and isothermal 
equations of state. In both cases the numerical diffusion introduced 
by the solver is sufficient to quench the development of instabilities. 
Section|4]deals with high resolution simulations of the development 
of these instabilities. 



3.1 Analytical approximations 

The overall structure of the colliding wind binary is given in Fig.[T] 
The density map shows two shocks separating the free winds from 
the shocked winds. The shocked winds from both stars are sepa- 
rated by a contact discontinuity. The Bernouilli relation is preserved 
across shocks hence 
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across the first shock. The subscript s refers to quantities in the 
shocked region and we have neglected the thermal pressure in the 
unshocked wind due to its high Mach number. A similar equa- 
tion holds for the second shock. The Bernouilli relation is con- 
stant in each shocked region but discontinuous at the CD. There, 
Pis = P2S by definition and v\ s = V2 S = on the line- 
of-centres so that the two Bernouilli equations combine to give 
pisvloi — p2 S v^ 2i with ps the value of the density on each side of 
the contact discontinuity. Assuming that the density is constant in 



each shocked region on the binary axis (the numerical simulations 
carried out below show this is a very good approximation) then 
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where pi (p2) is the value of the density at the first (second) shock. 
The above relation states the balance of ram pressures ([Stevens] 
|etal.|1992| >. Using Eqs. {]} and (6) then yields 
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where n is the distance between the first star and the first shock 
and r2, the distance between the second star and the second shock. 
If the shock is thin then n + r% ~ a and the distance R s « T2 of 
the CD to the second star is 
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Note that, for a given 77 ^ 1, the contact discontinuity is closer to 
the second star for a 2D geometry than for a 3D geometry. 

The shock positions are not easily derived away from the line- 
of-centres, where the density is not constant in the shocked winds. 
Analytic solutions have been derived based on the thin shell hy- 
pothesis, which considers both shocks and the contact discontinu- 
ity are merged into one single layer. |Stevens et al.](|1992|> (s ee also 
|Luo et al.|1990l|Dyson et al.|1993| and |^ntokfain et al.|2004) derive 
the following equation for the shape of the interaction region by 
assuming that it is located where the ram pressures normal to the 
shell balance: 
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The same analysis for the 2D structure (Eq.[6) leads to 
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|Canto et ah] ( |1996| ), extending the work of Wilkin] ( |1996| ), found 
an analytical solution in the thin shell limit based on momentum 
conservation (hence, taking into account the centrifugal correction 
i.e. the forces exerted on the gas as it follows a non-linear path 
along the shock, [Baranov et al.|197Tj|Dyson|1975l|Girard & Will-| 
|son|1987| >: 
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(see Fig.[T]for the definition of 0i and 02.) The same analysis in 2D 
leads to 
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3.2 2D study 

We performed a systematic study of the 2D geometry of the interac- 
tion zone in the adiabatic case for 77 ranging from 1 down to 1/128 
with Mach number M = 30 for both winds. Fig. [2] shows how 
the main features of the colliding wind binary vary with 77. The 
positions of the discontinuities on the binary axis (top left) were 
computed by determining the local extrema of the slope of the den- 
sity, excluding the masks. There is very good agreement with the 
analytic solution for the position of the standoff point (Eq.[l0|. The 
relation for the ratio of shock positions (Eq. [9} is also verified (top 
right). As 77 decreases both shocks and the contact discontinuity 
get closer to the star with the weaker wind. Since the thickness of 
the shell decreases as 77 decreases, proper modelling of the inter- 
action region for low 77 requires a higher numerical resolution. For 
77 < 0.25, the second wind is totally confined and there is a recon- 
finement shock on the line of centres behind the second star (see 
Fig. [TV This shock draws closer to the second star as 77 decreases 
(Fig^[2] bottom left). Similar structures were found by Myasnikov 
|& Zhekov| ( |1993| ) and |Bogovalov et al.|f2008| ) (in the latter case 
for 77 < 1/800). The last panel (bottom right) shows the asymp- 
totic opening angle of the contact discontinuity. The solution from 
|Stevens et al.| ( |1992] ) gives a better agreement for low values of 77. 

For given Mach numbers, the geometrical structure of the col- 
liding wind binary is set by 77. We performed a series of tests for 
77 = 1/8 = 0.125 and different combinations for Mi 
and M2 • Although the density and velocity fields were different in 
all cases, both shocks and the contact discontinuity were located 
at the same place along the line of centers. Further away from the 
star we notice that the reconfinement shock position changes up 
to ~ 25% when changing the velocity and mass loss rate of the 
winds. All other discontinuities are located at the same place. Sim- 
ulations with Mi — M2 — 100 do not show differences from the 
case Mi = M2 = 30, as could be expected since thermal pres- 
sure is negligible in both cases. However, the structure for given 
77 depends somewhat on the Mach number of the winds if these 
are not very large. Fig. [3] shows the density maps for 2D simula- 
tions with 77 = 0.25 but with different values for the wind Mach 
numbers obtained by changing the wind temperature. If both winds 
have M = 5 instead of M = 30, the shocked region is wider and 
a reconfinement shock appears at « 15a (beyond the region shown 
in Fig. [3}. The position of the contact discontinuity remains the 




Figure 2. Dependence of the shock geometry with 77 in 2D. Top left panel: 
Position of the different density jumps: first shock (black crosses), contact 
discontinuity (blue diamonds) and second shock (green asterisks). The 2D 
analytic solution for the contact discontinuity (Eq.[To) is overplotted (blue 
solid line). Top right panel: ratio of the shock positions measured from the 
simulations and compared to Eq.[9] Bottom left panel : position of the re- 
confinement shock. Bottom right panel: asymptotic opening angle (crosses) 
compared with the asymptotic angle derived from the Canto et al. (1996 1 
(dashed line) and Stevens et al. (1992 1 (solid line) solutions. 



same. When both winds have different Mach numbers, the whole 
shocked structure is more bent towards the wind with the higher 
Mach number : thermal pressure from the low Mach number wind 
is not negligible in the shock jump conditions (see Eq. [7} and the 
added term displaces the shock away from the low Mach number 
wind. 

We also investigated the overall structure in the isothermal 
case, quenching the strong instabilities that are present in this case 
(see §4.2| ) by using a highly diffusive solver. In this case pressure 
support is weaker and the shell is much thinner, as expected. The 
double shock structure and CD are only visible on the line of cen- 
tres when using a very high spatial resolution. The position of the 
thin shock structure on the line of centres is within 10% of the 
CD position found in adiabatic simulations. The asymptotic an- 
gle is difficult to assess as the shock structure is smoother than in 
the adiabatic case (see e.g. Fig. [4] below) but the bracketing val- 
ues are consistent with those found in the adiabatic case. We find 
that the weaker wind can be fully confined as in the adiabatic case. 
However, this occurs further away from the star than in the adia- 
batic case shown in Fig. [2] (at w 6.4a for 77 = 1/16 and 2.2a for 
77 = 1/32). 



3.3 3D study 

We completed this 2D study with the analysis of a few large scale 
3D simulations, computationally more expensive than the previous 
2D simulations. Fig. [4] shows the density maps for adiabatic and 
isothermal 3D simulations with 77 = 0.5 and 77 = 1/32 (,M=30). 
In the adiabatic case, one can clearly see the two shocks and the 
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Figure 3. Density maps for 2D simulations with 77 = 0.25 and different Mach numbers for the winds (.Mi, M.2). The 2D analytic solutions derived from 
the assumptions of Canto et al. (1996 1 and Stevens et al. (1992 1 are represented respectively by the dashed and solid line. The analytic solutions both assume 
infinite Mach numbers for both winds. 



contact discontinuity. For 77 = 1/32 the weaker wind is totally con- 
fined with maximum extension along the axis up to 5a away behind 
the star. For 77 = 1/64 « 0.016 (not shown) we find the reconfine- 
ment shock occurs at 1.0a. This is consistent with the 2D results 
(Fig. [2} if assuming the rough mapping ^Jfj^B —> ?72D suggested 
by Eq.|10| Indeed, we find no reconfinement shock for 3D simu- 
lations with 77 = 0.08 ( which would correspond to 772D ~ 0.29 
in Fig. [2]) Pittard & Dougherty (2006 ) performed 2D axisymmetric 
simulations showing a reconfinement shock for 
eta = 0.02 but not for 77 = 0.036. We performed several 3D sim- 
ulations with 77 = 1/32 = 0.03125 or 77 = 0.02 and for different 
values of the Mach number M (assumed identical in both winds). 
We found that reconfinement occurs in all cases when M = 30 or 
100 but that no reconfinement occurs for 77 = 0.02 or 77 = 1/32 
when M — 5. As in the 2D case, non-negligible thermal pressure 
has an impact on the structure of the colliding wind binary. Whereas 
the presence of reconfinement for low 77 and high Mach numbers 
around a threshold value 0.02-0.03 appears robust, the precise de- 
termination of this threshold value or of the properties of the re- 
confinement region is sensitive to the exact wind properties (Mach 
number). Radiative cooling, which is neglected here, can also have 
an impact on reconfinement (e.g. 2D isothermal simulation showed 
reconfinement further away from the star than in the adiabatic case, 
§3.2). 

The positions of the discontinuities along the line of centres 
agree within 2% with the expected values. As with the 2D case, 
the shock shape is better approximated by the solution of Stevens 
|et al.| ( [l992t at low 77. For 77 = 0.5 we find #oo = 71° whereas the 
asymptotic angle from both Steven s""et al.| ( [1992| and |Canto et al.| 
11996] ) give 78°; for 77 = 1/3 2 0.03125 we ge t 23° co mpared 
to theoretical estimates of 27° ( Stevens et al. 1992 ) and 35° ( |Canto| 
|et al.|[T99 6 ]>. On the ot her hand, Figs. 3j4 show that the analytic 
solution of Ca nto et aT] ( |1996| ) is a better approximation to the con- 
tact discontinuity shape at high 77. For 77 ^ 1/32, close to the line 
of centres, the shocked region is thinner in the 3D case than in the 
2D. For smaller values of 77, the shocked zone is thicker in the 3D 
case. In all cases the contact discontinuity is further away from the 
second star in the 3D case than in the 2D case. 

We have studied the geometry of the interaction region in 2D 
and 3D. We conclude that analytic solutions give satisfactory agree- 
ment with the results of the simulations. The solution based on ram 
pressure balance normal to the shock reproduces better the asymp- 
totic opening angle of the flow at low 77. We also find that the 



weaker wind can be entirely confined for low values of 77. How- 
ever, the interaction region is susceptible to instabilities that can 
modify these conclusions. This is investigated in the next section. 



4 INSTABILITIES 

4.1 The Kelvin-Helmholtz instability (KHI) 

When the exact Riemann solver is used, there is less numerical dif- 
fusion and the velocity shear at the contact discontinuity leads to 
the development of the KHI. The interface of two fluids is unsta- 
ble to any velocity perturbation along the flow in the absence of 
surface tension or gravity ( Chandrasekhar 1961 ). The growth rate 
of the instability in the linear phase is tkhi = A/ (2tvAv) where 
Av is the difference of velocity between the two layers and A the 
wavelength of the perturbation. In practice, numerical simulations 
are limited by diffusivity and the minimum resolvable structure, 
inevitably stunting the instability at small A. At the other end of 
the scale, the development of instabilities with large wavelengths 
can be hampered by their advection in the flow. The dynamical 
timescale can be estimated by Td yn ~ a/c s where c s is the post- 
shock sound speed, which is of the order of the wind velocity Voo 
in a strong adiabatic shock. Hence, the scale of the perturbations 
may be expected to be limited to A/a < Av/v. For two identical 
winds with terminal velocities of 2000 km s _1 and a — 1 AU, 
r d yn ^ 6.8 x 10 4 s= 2.2 x 10" 3 yr. 

We performed a set of simulations with 77 = 1, increasing the 
velocity Vooi of the first wind to investigate the impact of the KHI 
in the adiabatic case. The mass loss rate Mi was simultaneously de- 
creased and the Mach number 1 of the wind was kept equal to 30. 
The size of the domain is 8a and the resolution is n x — 128 with 5 
levels of refinement. The simulations were run up to t — 600rd yn - 
A steady state is reached well before the end of the simulation, as 
determined by looking at the time evolution of the total r.m.s. of 
the density or velocity perturbations over the whole simulation do- 
main. Restricting ourselves to this steady state interval, which we 
checked to be much longer than the advection time along the con- 
tact discontinuity, we then computed the time average of the veloc- 
ity r.m.s. for each cell of the domain. We used the median value 
over the same time period as our reference. The purpose was to 
quantify the saturation amplitude of the perturbations. 

The results are shown in Fig. [5] The upper panels gives the 
density maps for the different cases while the corresponding lower 
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Figure 4. Density maps for 3D simulations with 77 = 0.5 and 77 = 1/32 = 0. 
are located at the intersections of the dotted lines. The dashed line represents 
|et al. 1(1992) . The length scale is the binary separation a. 

panels show the time average of the r.m.s of the velocity fluctua- 
tions. No instabilities are present when the two winds are exactly 
identical, as expected since there is no velocity shear. Introducing 
a 10% difference in the velocity of the winds leads to low am- 
plitude perturbations that are significant only close to the contact 
discontinuity. A dominant wavelength can be identified, probably 
because growth for such a weak velocity shear is restricted to a 
small domain by diffusivity at short wavelengths and advection at 
long wavelengths. The r.m.s. of the velocity and density pertur- 
bations saturates at about 10%. When small scale 
eddies are visible. They are stretched in the direction of the flow. 
The position of the shocks is barely affected by the instability. The 
perturbations affect a larger zone on both sides of the contact dis- 
continuity but their amplitude remains around a few tens of per- 
cent r.m.s., somewhat higher for the density than for the velocity 
perturbations. When Vooi — 20voo2 (fourth panel) the instability 
has become non-linear judging by the 100% r.m.s. of the velocity 
(and density) fluctuations. The location of the contact discontinu- 
ity fluctuates significantly yet the region with the strongest r.m.s. is 
not much wider than for the previous cases. We also investigated 
in this last case whether keeping the wind temperature constant as 
Vooi is varied, instead of keeping Mi constant, led to differences. 
The outcome was similar. 

A similar set of simulations was performed with 77 = 1/16 = 
0.0625 (Fig [6). There is no velocity shear or contact discontinuity 
when Vooi = V002, even in the case rj ^ 1. This can be proven as 
follows. The Bernouilli constant (Eq.[7} has the same value in both 
shocked region when Vooi = V002, so the densities are identical at 
the contact discontinuity (where pressures equalise) on the line-of- 
centres. The gas is polytropic with P = Kp~ J and K constant in 
each region. Writing that p and P are equal on both sides of the 
contact discontinuity on the line-of-centres requires that K has the 
same value in both shocked regions. Therefore, p± s = p2s along 
the contact discontinuity. Using that the Bernouilli constant is the 
same in both shocked regions then proves that vi s = V2 S at the con- 
tact discontinuity. Actually, there is no discontinuity in this case. 
The simulation with Vooi = V002 confirms that there is no velocity 
shear and that the KHI does not develop. When Vooi = I.IV002 
only weak perturbations are seen, limited to a small region close 
to the contact discontinuity. A dominant wavelength can be iden- 
tified as in the case rj = 1. When Vooi = 2voo2 the center line 
of the perturbations approximately matches the shape of the unper- 
turbed contact discontinuity. The first shock is not affected by the 
instability. The velocity perturbations affect all the region of the 
shocked second wind and part of the shocked wind of the first star. 
The density perturbations have a higher r.m.s. than the velocity per- 
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turbations, reaching close to 100% close to the contact discontinu- 
ity. The velocity perturbation are strong when v^i = 20^002 and 
are mostly confined to the shocked second wind. High r.m.s. den- 
sity fluctuations extend to the first wind, distorting slightly the first 
shock. (The sawtooth appearance of the wings in the vioo — V20 00 
r.m.s. maps are an artefact of the limited time range over which the 
average was done.) The backward reconfinement of the wind of the 
second star is affected by the instability, occurring much closer to 
the second star than in the case with equal wind velocities. 

The KHI modifies the interaction region as soon as the wind 
velocities are slightly different. The simulations suggest that the 
relative amplitude of the perturbations becomes significant when 
vi > 2^2, although we cannot rule out that limited numerical res- 
olution does not impact the growth of the instability for smaller 
velocity shears. The instability does not erase completely the con- 
tact discontinuity. However, the turbulent motions tend to smooth 
out the initial structures in the region of the wind with the smaller 
velocity. 



4.2 Isothermal equation of state: thin shell instabilities 

When thermal support in the shocked zone is too weak, the shell 
becomes thin and unstable. This occurs for instance when the adia- 
batic index is decreased ( Mac Low & Norman 1993). More realis- 
tic numerical simulations including radiative cooling functions also 
show the shocks become thinner and unstable as cooling increases 
( [Stevens et al.|1992||Pittard|2009| but see |Myasnikov et al.|199"8] >. 
The instability is usually referred to as 'thin shell instability' al- 
though several physical mechanisms may be at work, including the 
KHI. The non-linear thin shell instability (NTSI, | Vishniac| 1 994| ) is 
found in hydrodynamical simulations when the thin shell is moved 
away from its rest positions by perturbations with an amplitude at 
least greater than the shell width (Blondin & Marks 1996). The in- 
stability is due to an imbalance in the momentum flux within the 
shell as shocked fluid moves towards opposing kinks. The trans- 
verse acceleration instability (TAI, Dgani et al. 1993 1996) occurs 
when at least one of the colliding flows is divergent and assumes an 
infinitely thin shell. Both linearly unstable breathing and bending 
modes are found. The breathing mode is due to the acceleration of 
the flow along the shell whereas the bending mode arises from the 
mismatch in ram pressure of the wind impacting each side of the 
thin shell when it is displaced from its equilibrium value. 

We studied the growth of thin shell instabilities in colliding 
wind binaries using 2D simulations with an isothermal equation 
of state. Initial investigations showed that the thin shock structure 
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Figure 6. Same as Fig.[5]but for 77 = 1/16 = 0.0625. 
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( §3.2| ) becomes unstable only if there are a sufficient number of 
cells available (> 4) to resolve the shock structure. The minimum 
number of cells required is even larger if a highly diffusive solver is 
used. Low resolution simulations without mesh refinement (256 x 
256 cells) do not resolve the shock structure and stay stable. We 
decided to use those steady state solutions as the initial input for 
simulations at higher resolution, so as to be able to study in as much 
as possible the initial linear growth phase of the instabilities. The 
winds are chosen to have identical velocities in order to exclude 
any seeding by the KHI (jp~T). 

The evolution of a colliding wind binary with rj = 1, identi- 
cal velocities and an isothermal equation of state is shown in Fig. [7] 
The size of the domain is 3a. The left panels show the case with one 
level of mesh refinement, the right panels show the case with four 
levels. At low resolution (left panels), perturbations become visible 
away from the line-of-centres early in the simulation (t = 9.5 x 10 4 
s). These perturbations grow slowly as they are advected, thicken- 
ing the layer. At t = 1.5 x 10 5 s another instability develops close 
to the binary with a growth rate faster than the advection rate and 
a distinct morphology. In this case matter piles up in the convex 
parts of the shell, which move steadily away from the initial shock 
position without the oscillatory behaviour seen in the wings. At the 
end of the simulation it — 3.1 x 10 5 s) the colliding wind region 
is dominated by these large scale perturbations. At higher resolu- 
tion (right panels), the initial instability appears earlier and is also 
present closer to the binary axis. At £ = 9.5 x 10 4 s there already 
is a superposition of modes and one cannot define a unique wave- 
length any more. At t = 1.8 x 10 5 s oscillations are present even 
on the binary axis and the structure is not symmetric any more. The 
final density maps shows a thicker shell with small scale structures. 
The oscillations are smaller than for the low resolution simulation 
at this time. The evolution at subsequent times shows comparable 
amplitudes in the oscillations at high and low resolution. 

Similar behaviour was described by Blo ndin & Koerwer] 
(1998 ) in their simulations of stellar wind bow shocks. We ten- 
tatively associate the small amplitude instability that develops first, 
away from the binary axis, with the TAI. This is a linear instability 
that can be seeded by the initial numerical noise. The large ampli- 
tude instability that develops later on the binary axis is likely to be 
the NTSI. We examine below the supporting evidence. 

4.2. 1 Evidence for the Non-linear Thin Shell Instability (NTSI) 

The NTSI shows the highest growth rate for perturbations of order 
of the shell wi dth L. The the oretical estimate is r t h = L/c s = 
2.0 x 10 4 s (Vishniac 1994) for the parameters appropriate to 
our simulations, smaller than the advection timescale (Td yn — 
6.8 x 10 4 , increasing near the binary axis as the flow velocity in 
the shocked region goes to zero on axis). Hence, the fastest growing 
mode of the NTSI should be seen, independently of the numerical 
resolution, as long as the shell is resolved. We compared this esti- 
mate with the time evolution of the velocity perturbations in four 
simulations with 1,2,3 and 4 levels of refinement, using an exact 
Riemann solver. For each simulation, we computed the r.m.s. of the 
velocity for a line of cells along the binary axis, where the NTSI is 
presumed to dominate. We normalised the data to the value at the 
same arbitrary reference time taken close to the beginning of each 
simulation. The r.m.s. were smoothed to suppress small wavenum- 
ber perturbations that appear at high resolutions. The logarithm of 
the r.m.s is shown in Fig. [8] The shell readjusts to the higher nu- 
merical resolution up to t ~ 9.5 x 10 4 s. Close inspection of the 
density maps reveals the presence of density fluctuations on the 
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Figure 7. Density maps showing the evolution of a 2D colliding wing binary 
when 77 = 1 and 7 = 1.01. Time is given in seconds, t = corresponds 
to the restart at high resolution of an initial low resolution simulation (256 
cells, no mesh refinement). On the left panel there is one level of refinement 
(maximum resolution equivalent to 512 cells), on the right panel there are 
four levels (maximum resolution equivalent to 4096 cells). 



scale of the shock width during this transition. This numerical re- 
laxation triggers the NTSI close to the binary axis (left panels of 
Fig. [7]). In the simulations with highest numerical resolution (right 
panels of Fig. [7} the NTSI develops in regions that are already per- 
turbed by the growth of the first instability (most likely the TAI, see 
j ]4.2.2| l. These fast growing perturbations may contribute to trigger 
the NTSI. The NTSI moves the shock away from its rest position as 
the bending modes are amplified and mass collects at the extrema 



Unstable colliding winds at high resolution 9 




Figure 8. Logarithm of the root mean square of the velocity on the line of 
centres as a function of time. The curves represent maximal resolutions of 
512 (dotted), 1024 (dot-dashed), 2048 (dashed) and 4096 (solid) cells per 
dimension. The thin straight lines show the fits to the linear phase for each 
resolution. 



very short. Another is that we found that our initial velocity pro- 
file along the shock is inconsistent with the equilibrium solution 
proposed by Dgani et al. (1993). This was corrected by Myasnikov 



et al. (|1998| > but they concluded that the set of equations used by 
Dgani et al. (1993) led to inconsistencies in the dispersion relations, 



casting doubt on the theoretical rates to expect. We suggest that it 
is not possible to neglect, as was done, the derivatives d/dO in the 
equations (0 corresponds to the polar angle to the binary axis with 
the origin at the stagnation point), since there is a significant change 
in the azimuthal speed of the incoming flow as it is decelerated and 
redirected along the shock. Although our results still support the 
presence in the simulations of some form of the TAI, the simu- 
lations also show that the saturation amplitude of this instability is 
low compared to the NTSI. In all the simulations we performed, the 
non-linear evolution was dominated by the large scale, high ampli- 
tude perturbations induced by the NTSI. At best, the TAI may play 
a role in the early stages as a seed instability for the NTSI, as de- 
scribed in §4~2A\ 



( Vishniac 1994). The exponential growth timescale estimated from 
fitting the r.m.s values are r « 3.1 x 10 4 , 2.9xl0 4 , 4.5xl0 4 and 
4.7 xlO 4 s for increasing resolutions (mesh refinement). There is 
an increase of 50% of the measured growth timescale whereas the 
cell size (and therefore the available wavelength range potentially 
accessible) increases by a factor 16. This is in reasonable agreement 
with the theoretical value and the expected behaviour with chang- 
ing resolution, confirming that the NTSI is triggered in our simula- 
tions. Fig. [8] also shows that the saturation amplitude is somewhat 
smaller as the resolution is increased (compare also the bottom left 
and right panels of Fig. [7} and that it converges to a resolution- 
independent value. 



4.2.2 Evidence for the Transverse Acceleration Instability (TAI) 

The numerical simulations show that the initial perturbations are 
preferentially located off the binary axis, have an oscillatory be- 
haviour with a small wavelength and grow faster when the spatial 
resolution is increased (Fig. [7}. The rapid development of these 
perturbations is consistent with a linear instability. These proper- 
ties are reminiscent of the TAI. The TAI studied by Dgani et al.| 
( 1993 1996) is an overstability with an oscillation frequency of the 
velocity perturbations oc 1/A. The growth timescale is oc y/X and 
indeed smaller wavelength perturbations grow faster at higher reso- 
lution. Vishniac (1994) noted that the growth is limited by pressure 
effects and that the TAI grows faster than the NTSI when 



Here, I is the minimum distance along the contact discontinuity 
(I = on the binary axis) beyond which the TAI can develop for 
a given wavelength A. The relevant wavelengths are smaller than 
R s and larger than the shell width L ~ R s /M 2 , with the smaller 
scales growing faster. The instability develops preferentially along 
the wings jBlondin & Koer wer 1998). The presence of the TAI 
closer to the binary axis at the highest resolution may explain why 
the growth rate of the NTSI (see Fig. [8) does not perfectly match 
the theoretical value. 

Despite the similarities, we could not formally identify the 
TAI. One difficulty is that we were not able to quantify the growth 
rates as several modes interact quickly and make the linear phase 



4.2.3 Evolution with an initial velocity shear and at low r\ 

In real systems the velocities of the winds are never exactly equal 
and the contact discontinuity is subject to the KHI. Even for a 1% 
velocity difference between the winds, this instability theoretically 
has a larger growth rate than the TAI and NTSI. Fig. [9] compares 
simulations for r\ — 1 with equal winds or v\oo — 2v<2oo, subject 
to the KHI. We also include here a map of the r.m.s. of the velocity 
fluctuations observed over a long averaging period. There is little 
difference in the outcome between equal winds and v\oo — 2v2oo, 
either in the appearance of the turbulent region (top row) or in the 
r.m.s. of the perturbations (second row). If anything, the KHI seems 
to increase slightly the region where strong fluctuations occur. The 
NTSI dominates the final non-linear phase even when the KHI is 
initially present. The r.m.s. values close to one are the expected 
outcome of the NTSI ( [Vishniac|1994] ). 

We found the same results for simulations with 77 = 1/16 = 
0.0625. The corresponding density maps and velocity perturbations 
are given in the bottom two rows of Fig. [9] The NTSI was stud- 
ied theoretically for planar shocks but the 77 = 0.0625 simulations 
show it is also present and dominant when the shock is curved, 
although following it requires high numerical resolutions. The sim- 
ulations were performed with n x = 128 and 5 levels of refinement 
in a box of size 8a. For lower resolutions the NTSI is not triggered 
and the final result is stable (the same is observed for 77 = 1). The 
density maps for equal winds and vioo = 2v 2o o look similar. The 
highest velocity perturbations are at the same location but the r.m.s 
values are higher when an initial shear is present. We conclude that 
having a velocity shear in a thin shell increases the amplitude of the 
perturbations but does not affect much the morphology of the un- 
stable flow, which is mostly set by the NTSI. This is consistent with 
Blondin & Marks (1996) who concluded from their simulations of 
perturbed slabs that the KHI does not strongly modify the outcome 
of the NTSI. 



4.2.4 Effect of increasing pressure in the stellar winds 

Pressure has a stabilising effect on both instabilities. We performed 
a simulation with M 1 =M 2 =6 with all other physical and numeri- 
cal parameters identical to those of the 77 = 1, vioo = V200 simula- 
tions. Both instabilities are seen to develop but more slowly. Keep- 
ing the wind velocity constant, a lower Mach number implies a 
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Figure 9. Top row: density maps for 77 = 1 with tj 1oo =- tj 2 oo (left panel, 
from the same model shown in Fig. [7} and vioo = 2tj2oo (right panel). 
Second row: corresponding time-averaged r.m.s. of the velocity fluctuations 
(on a log scale). Bottom two rows: same for 77 = 1/16 = 0.0625. 



higher sound speed but the thickness of the shell increases faster 
so that the growth timescale of the NTSI (oc L/c s oc l/M) is 
longer. The NTSI is also harder to trigger as it requires a per- 
turbation of amplitude comparable to the size of the shell. The 
TAI develops more slowly as pressure suppresses the development 
of small wavelengths perturbations in the radial directions ( Dgani 
|et al.|1993) . The final non-linear phase with high amplitude pertur- 
bations, shown in Fig. [To] appears later than in Fig. [7] The shell is 
indeed thicker and presents smaller density contrasts than for high 
Mach numbers. Comparing Fig. [9] with Fig. [To] the amplitude of 
the variations in shock location or the r.m.s of the fluctuations do 
not appear to change much but the oscillations in shock location 
seem to have a longer wavelength. 



Figure 10. Left: density map of 2D colliding wing binary when 77 = 1 , 
7 = 1.01 and A4 = 6 for the highest resolution. Time is given in seconds. 
Right: time-averaged map of the r.m.s of the velocity fluctuations. 



4.3 A comparison of unstable adiabatic and isothermal cases 

Finally, we compare the non-linear outcome of simulations with 
unstable colliding wind regions in the isothermal and adiabatic 
cases. Figs. 5j6 and Fig. [9] show cases with 77 = lor?7 = 1/16 
and vioo = 2^2oo for both the adiabatic and isothermal cases. The 
r.m.s. amplitude is larger for isothermal winds than for adiabatic 
winds when the same wind parameters are used. The unstable re- 
gion extends beyond the wings of the contact discontinuity in the 
case of isothermal winds, unlike the adiabatic case where most of 
the fluctuations seem to take place within the shocked region of 
the weaker wind. The NTSI creates more small scale structures and 
higher density contrasts are possible when the winds are isother- 
mal. The weaker wind still propagates freely over a significant frac- 
tion of the domain despite the strong perturbations at the interface 
in the isothermal case. In contrast, the adiabatic simulations show 
that the free flowing weaker wind is confined to a very small region 
(j ]4.1| ». The wind is still expected to be confined at some distance 
from the star in the isothermal case (see §3.2) but this happens fur- 
ther away than in the adiabatic case even when the thin shell insta- 
bilities develop. 



5 DISCUSSION 

5.1 Morphology of the interaction region 

We have carried out 2D and 3D hydrodynamical simulations of col- 
liding winds to study the morphology of the interaction region and 
the instabilities that can affect it when orbital motion can be ne- 
glected. We first examined the relevance of widely-used analytical 
estimates. The position of the standoff point is very well predicted 
by the standard ram pressure balance on the line-of-centres. Away 
from the binary axis, when 77 is close to 1, the opening angle of the 
contact discontinuity is well approximated by the analytical solu- 
tion proposed by |Canto et al.|(T996] >, which assumes conservation 
of mass and momentum in a thin shell. The semi- analytical solu- 
tion of Steven s et al.| ( |1992| >, which assumes balance of the ram 
pressures normal to the surface, is a better approximation when 
77 <C 1. This clarifies the range of validity for these approximations 
that have found widespread practical use in the literature. 

Numerical simulations also show that the weaker wind can 
be fully confined for low 77, with the presence of a backward ter- 
mination (reconfinement) shock, for both isothermal and adiabatic 
winds. The region where the weaker wind propagates freely is re- 
duced when the Mach number of the wind is small, when the KHI 
develops or when the wind is isothermal. This may have some ob- 
servational consequences. One possibility is that the lines from the 
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confined wind show unusual profiles or intensities because the wind 
terminates very close to the star. Another possibility is stronger, 
variable absorption instead of smooth absorption when the line-of- 
sight crosses the region where a freely-expanding wind is expected. 

More realistic simulations would include wind acceleration 
and radiative inhibition or brak ing ([Stevens & Pollock [1994 



|Owocki & Gayley|[T995l |Pittard||2009[ |Parkin & Gosset 



2011) 



The wind velocity at the stagnation point is then different from 
its asymptotic value, increasingly so when ram pressure balance 
occurs close to one of the stars. The principal consequence is to 
change the location of the stagnation point ( Antokhin et al. 2004). 
The basic geometry of the interaction region does not change al- 
though the asymptotic values e.g. of the contact discontinuity are 
probably best described by some effective r\. In some extreme cases 
a stable balance may not be achieved and the wind- wind collision 
region collapses onto the star with the weaker wind ( Stevens et al. 
[T992) |PittardTT998] > . Another possible consequence is that a veloc- 
ity shear may appear even if the coasting velocities of the winds 
are assumed to be equal, generating the KHI where it would not be 
expected. 

Orbital motion must be included when studying the large-scale 
structure of colliding winds. The interaction region wraps around 
the binary at distances of order Vooftrb, where Voo is the veloc- 
ity of the stronger wind (Walder et al. 1999). On smaller scales 
(intra-binary), a non-zero orbital velocity skews the interaction re- 
gion by an angle tan a ~ v rb / Voo at the apex ([Parkin & Pittard 



2008). The opening angles of the shocks are slightly modified on 
the leading and trailing edges but the morphology of the interaction 
region does not dramatically change on scales <C ^oo^orb (Lemas- 
|ter et al.|[2 007). Exploratory simulations show that the reconfine- 
ment shock is still present when orbital motion is included in a low 
77 model. According to our results ( §3.3| >, no such shock is expected 
to form in the adiabatic simulation of van Marie et al. ( 2011 ) since 
it has 77 = 1/7.5 w 0.14. Reconfinement shocks can occur at some 
phases and not at others in binaries with highly eccentric orbits, as 
different cooling or wind velocities are probed when the separation 
changes (e.g. the periastron passage of the 77 Carinae, see Parkin 
|et al.|201l) . The morphology also depends on the history of the 
shocked gas and can exhibit strong hysteresis effects in eccentric 
systems ( |Pittard|2009] >. 



5.2 Impact of instabilities 

Hydrodynamical instabilities have a major impact on the structure 
of the colliding wind binary. Although the overall aspect of the in- 
teraction region can still be recognised in a time-averaged sense, the 
wind interface can become highly turbulent, generating strong time 
and location-dependent fluctuations in the flow quantities. Velocity 
shear at the contact discontinuity in the shock region leads to the 
development of the KHI. An accurate Riemann solver is required 
to follow this instability. Eddies are already present at the interface 
even with a 10% velocity difference. The amplitude of the perturba- 
tions can be significant with r.m.s. values in the tens of percent for 
the case of adiabatic colliding winds with vioo = 2v2oo. The mix- 
ing is limited to the region of the weaker wind, with the strongest 
perturbations located close to the initial contact discontinuity. The 
KHI has no impact on the location of the stagnation point. Equal 
winds are not expected to trigger the instability but introducing or- 
bital motion was found to generate a small velocity shear even for 
this case ( [Lemaster et al.| 2007 ). Curiously, | van Marie et al.| ( [20TT| ) 
find the opposite i.e. no KHI for nearly adiabatic winds with orbital 
motion, v\oo = 1.3^2oo and 77 = 0.14. We would expect to see 



significant mixing in the inner binary system, where the interaction 
region is only slightly skewed, unless it is dampened by numerical 
diffusion. 

In isothermal simulations, an instability reminiscent of the 
TAI develops initially away from the binary axis. A second insta- 
bility develops on the axis whose growth rate and properties iden- 
tify as the NTSI. The NTSI dominates the non-linear evolution of 
isothermal colliding winds, leading to highly turbulent structures 
and large amplitude fluctuations in the location of the interface, 
including the stagnation point on the binary axis. Our results con- 
firm the conclusions of Blondin & Koerwer (1998) who stressed 
the dominance of the NTSI and the stabilising effect of pressure 
in their simulations of bow shocks. They also saw 'wiggles' devel- 
oping early on in the shock with the same properties as those we 
attribute to the TAI-like instability. The trigger for the NTSI is not 
discussed but it is likely provided by the wiggles. However, they 
did not attribute these to the TAI and instead argued that the TAI 
acts only once the shell is perturbed by the NTSI. 

The presence of instabilities in real systems is probably un- 
avoidable. The KHI may lead to moderate mixing of the material in 
adiabatic situations. The strongest mixing is obtained for high ve- 
locity shears which, in astrophysical systems, is likely to mean that 
at least one of the winds is radiatively efficient and not adiabatic. 
The radiative efficiency of the wind is classically parametrized by 
the ratio \ of the cooling and advection timescales, which can be 
evaluated as (Stevens et al.|1992] > 
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with x ^ 3 for an adiabatic wind and % < 3 for a radiatively ef- 
ficient wind. The ratio X1/X2 is therefore oc (vi/v2) 5 rj. Because 
v appears with a large power, a significant difference in wind ve- 
locities essentially implies that the slowest wind will be close to 
isothermal. In this case, thin shell instabilities develop but their 
outcome may be different because of the stabilising effect of ther- 
mal pressure from the neighbouring adiaba tic shock ([Stevens et ah] 
[T992l[WaTder & Folini| 1998) |Pittard|2009l parkin & Pittard|2010| 
van Marle^taLpOlT}. For thin, highly radiative shocks, the NTSI 
can probably be triggered by wind variability or changes in shock 
width as x varies along the orbit, if it is not already triggered by 
the TAI or KHI. The saturation amplitude depends strongly on the 
radiative losses and including a realistic cooling function in the en- 
ergy equation of the fluid is essential for a detailed comparison with 
observations ( [Strickland & Blondin| 1995) | Walder & Folini|1996) . 
The shock will necessarily be larger than the idealised isothermal 
case so the saturation amplitudes of the fluctuations can be expected 
to be in between the adiabatic and isothermal cases. Other instabil- 
ities may also be at work in radiative shells ( Chevalier & Imamura 
1982, Wal der & Folini|19 96). Compressed magnetic fields in the 
shock region, if present, can also modify the growth rates and satu- 
ration amplitudes. For instance, the KHI is stabilised when the flow 
is parallel to the magnetic field and the velocity shear is smaller 
than the Alven speed ( |Gerwin|1968] ). |Heitsch et al.| ( |2007) find that 
an ordered magnetic field has a stabilising effect on the NTSI in a 
thin slab. 

In conclusion, the impact of the instabilities studied here is 
conveniently summarised by saying that some amount of variability 
and mixing is expected in all cases but that the strongest variability 
and mixing are expected to be associated with the most radiative 
(hence luminous) colliding winds. 
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5.3 Computational requirements 

Following these instabilities is computationally demanding, espe- 
cially for low momentum flux ratios 77, and imposes a minimum 
spatial resolution together with an accurate Riemann solver. There 
are three numerical constraints on the spatial resolution. First, there 
must be enough cells within the stellar masks to properly generate 
the winds. For a coasting wind the mask can be larger than the ac- 
tual size of the star. This cannot be the case if the stagnation point 
is close to one of the stars low 77) and/or if wind acceleration, brak- 
ing or inhibition is taken into account. The second condition is that 
the resolution must be sufficient to resolve the location of the stag- 
nation point on the binary axis. This is increasingly demanding as 
77 decreases, but the increase in computational cost is steeper when 
working in the 2D setup (see §3.1| ). The last conditions relates di- 
rectly to the instabilities. For 77 = 1/32 = 0.03125, in a 8a simula- 
tion box, we found that a simulation with n x — 128 needs 7 levels 
of refinement in order to avoid numerical damping of the insta- 
bilities. At lower resolutions we see the initial development of the 
TAI far from the binary but it is quickly advected out of the sim- 
ulation box without being maintained. The NTSI is not triggered 
and the final result is stable. We find that the shell needs to be re- 
solved by at least 4 computational cells on the binary axis in order 
to develop the NTSI. Resolving the shell i.e. shock structure is the 
stringiest constraint on the numerical resolution. The thickness of 
the shell for the 2D adiabatic simulations given in Fig.[2](upper left 
panel) can be used to estimate the numerical resolution required to 
achieve this for a given 77. It drastically decreases for low values 
of 77 (slightly less so in 3D, which show thicker structures when 
77 <c 1/32 = 0.03125, see ^33). The shell width is thinner in the 
isothermal case so the values derived from Fig. [2] are strict lower 
limits for the required resolution. 

Large scale simulation of a system with low 77 and isother- 
mal winds require high resolutions for the instabilities to develop. 
The NTSI develops at slightly lower resolutions when the KHI is 
present and acts as the initial seed perturbation. For instance, with 
77 = 1/32, isothermal winds and v\oo = 2voo the NTSI develops 
with 6 levels of refinement instead of 7 in the case of equal winds. 
However, it seems that the effect decreases with lower values of 
77. The shell always needs to be resolved, if only minimally, be- 
cause the NTSI involves an imbalance of momentum within the thin 
shock layer. The Kelvin-Helmholtz instability in adiabatic winds is 
easier to model. It develops even for low resolution simulations 
when the velocity difference between both winds is large enough. 
For 77 = 1/32, adiabatic winds and the instability 

develops for 4 levels of refinement. The study of the large scale 3D 
evolution of unstable colliding winds remains a tremendous com- 
putational challenge. 



6 CONCLUSION 

We have studied the morphology and the instability of colliding 
wind regions using numerical simulations. Compared to previous 
works, our study extends to much lower values of the wind mo- 
mentum ratio, larger simulation domain and higher spatial resolu- 
tion thanks to adaptive mesh refinement. We investigate the appli- 
cability of semi- analytical estimates for the contact discontinuity, 
finding that the solution of Ste vens et al.| ( |1992) is the best approx- 
imation to the asymptotic opening angle for small 77. We find that 
the weaker wind can be entirely confined to a small region instead 
of expanding freely up to infinity over some solid angle when low 



77 colliding winds are considered in both the isothermal and adi- 
abatic limits. Instabilities in the colliding wind region are impor- 
tant because of the mixing and variability they induce. Resolving 
the shock structure is required to follow the development of insta- 
bilities, which imposes increasingly stringent minimal numerical 
requirements for smaller 77. Simulations that do not meet these re- 
quirements artificially dampen the instabilities that may be present. 
We follow the evolution of the KHI triggered by the velocity shear 
at the contact discontinuity between two winds and show that the 
eddies yield large fluctuations even for moderate initial shears. We 
formally identify the NTSI in our isothermal simulations and find 
that it dominates the long-term behaviour. Another instability, sim- 
ilar to the TAI, is present at the beginning of the simulations. Thin 
shell instabilities yield large fluctuations of the flow quantities over 
a wide region. Our study clarifies several issues in colliding wind 
binary models and provides a basic framework to which the results 
of more complex simulations, including additional physical effects, 
can be compared. 
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